Numerical analysis of backreaction in acoustic black holes 
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Using methods of Quantum Field Theory in curved spacetime, the first order in h quantum 
corrections to the motion of a fluid in an acoustic black hole configuration are numerically computed. 
These corrections arise from the non linear backreaction of the emitted phonons. Time dependent 
(isolated system) and equilibrium configurations (hole in a sonic cavity) are both analyzed. 



PACS numbers: 04.62.+V, 04.70.Dy, 47.40.Ki 
I. INTRODUCTION 

Black hole radiation predicted by Hawking in 1974 p] 
is one of the most spectacular results of modern theoret- 
ical physics. 

Even more surprising is the fact that this effect is not 
peculiar of gravitational physics, but is also expected in 
many completely different contexts of condensed matter 
physics HEI1- A fluid undergoing supersonic motion is 
the simplest example of what one calls an "acoustic black 
hole" . For this configuration Unruh , using Hawking 
arguments, predicted an emission of thermal phonons. 
This emission affects the behaviour of the underlying 
fluid because of the non linearity of the hydrodynami- 
cal equations governing its motion. 

Using methods borrowed from Quantum Field Theory 
in curved spacetime, this quantum backreaction has been 
studied for the first time in [6( , where the the first order in 
h corrections to the classical hydrodynamical equations 
were given. Because of intrinsic mathematical difficul- 
ties, the analysis was restricted to the region very close 
to the "sonic horizon" of the acoustic black hole; i.e., the 
region where the fluid motion changes from subsonic to 
supersonic. There, analytical expressions for the quan- 
tum corrections to the density and velocity of the mean 
flow have been provided. 

However, to have a detailed description throughout the 
entire system one has to proceed with numerics. This will 
be the aim of our present paper, which is organized as 
follows: in Sec. |H] we outline the classical fluid configu- 
ration which describes an acoustic black hole; the quan- 
tum backreaction equations are discussed in Sec. lIIII with 
emphasis on the choice of quantum state in which the 
phonons field has to be quantized. In Sees . II VI and IVl we 
give the numerical estimates for the quantum correction 
to the mean flow in two different cases: isolated system 
and system in equilibrium in a sonic cavity respectively. 
Section IVTl contains the final discussion. 



II. THE ACOUSTIC BLACK HOLE 

An acoustic black hole is a region of a fluid where its 
motion is supersonic. Here sound can not escape up- 
stream being dragged by the fluid. The boundary of this 
region is formed by sonic points where the speed of the 
fluid equals the local speed of sound. This is the acoustic 
horizon. A simple device to establish a transonic flow is 
a converging diverging de Laval nozzle 0,0. For station- 
ary free fluid flow the acoustic horizon occurs exactly at 
the waist of the nozzle. 

The basic equations describing the system at the clas- 
sical level are the continuity and the Bernoulli equations. 
We assume a one-dimensional stationary flow, therefore 
all the relevant quantities depend on z only, the spatial 
coordinate running along the axis of the de Laval nozzle. 
The continuity equation then reads 

A(z)p{z)v{z) = const = D , (1) 

where A is the area of the transverse section of the nozzle, 
p the fluid density and v the fluid velocity. The Bernoulli 
equation, under the above hypothesis, gives 

y+M/>) = 0, (2) 

where p(p) is the enthalpy. We have further assumed the 
fluid to be homentropic and irrotational. The speed of 
sound c is defined as 

c = p— . (3) 

dp 

For constant c (the case we consider) integration of © 
gives 

M (p)=c 2 ln^, (4) 
Po 

which inserted in Bernoulli equation yields 

v 2 

p = p Q e~^ , (5) 
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FIG. 1: Dependence of the cross section of the de Laval nozzle 
on the position z for a velocity field given by Eq. © and 
depicted in Fig. |5] The vertical dashed line corresponds to 
the location of the sonic horizon zh = 0.025 cm. 
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FIG. 2: Velocity (top) and density (bottom) from Eq. ® . The 
vertical dashed line correspond to the location of the sonic 
horizon (zh = 0.025 cm). The sonic black hole corresponds 
to z < zh- 



with po a constant. The assumed constancy of the speed 
of sound also gives the pressure p as p = c 2 p. The velocity 
profile describing the acoustic black hole is chosen as 



v = c < — arctan[/3(z — Zh)] — 1 

7T 



(6) 



where z = zh denotes the position of the waist of the 
nozzle (the sonic horizon). In the laboratory frame the 
fluid is moving from right to left, so v < and the sonic 
horizon occurs where v = — c. The constant D entering 
the continuity equation is determined by requiring the 
fluid to be sonic at the waist; i.e., 



D = -cA H p e 



-1/2 = _ 



ArPh 



(7) 



where Ah is the area at the horizon, pn the pressure 
and po — phg 1 ^ 2 / c 2 . Given this, the profile of the nozzle 
can be computed from Eq. (JIJ and is depicted in Fig. ^ 
where we have used Ah = 10~ 8 cm 2 , (3 — 600 cm -1 , 
Ph = 2 x 10 6 Pa and c = 250 m/s. These latter two are 
typical values for liquid Helium. The profiles of density 
and velocity are shown in Fig. [21 where the significant 
range of z is [0,0.05] cm; the horizon lies at zh = 0.025 
cm and its location is indicated by a vertical dashed line 
in the figures. In the region z > zh the motion of the 
fluid is subsonic; the acoustic black hole is the region 

z < ZH- 

As shown by Unruh (2| , sound waves propagating in an 
inhomogeneous fluid are described as a massless scalar 
field propagating in an effective curved spacetime de- 
scribed by an "acoustic metric" which depends on p and 



v. Quantization of these modes leads to the conclusion 
that in presence of a sonic horizon a thermal emission of 
phonons is expected, in complete analogy of what Hawk- 
ing found for gravitational black holes. The emission 
temperature of the phonons is Tjj = hk/ (2-kckb), where 
k is the surface gravity of the sonic horizon, defined as 



2 dn K 



(8) 



Kb is the Boltzmann constant and n is the normal to the 
horizon. For the specific acoustic black hole model we 
consider in this paper, T H = 1.1598 x 10~ 5 °K. 



III. THE BACKRECTION EQUATIONS 

The phonons quantum emission previously discussed 
modifies the underlying fluid flow according to the back- 
reaction equations derived in Ref. jf|, to which we refer 
for further details. For a one-dimensional flow, they read 



Ap B + d z {Ap B v B ) = d z 



A 



.( T (2)) . 



(9) 
(10) 



Here ps and vb are the quantum corrected density and 
velocity fields and ip B is the velocity potential; i.e., 
d z tpB — v B ; the overdot stands for time derivative. The 
(T ab ) which drive the backreaction are the quantum ex- 
pectation values of the pseudo energy momentum tensor 




FIG. 3: Quantum sources in the Unruh state (solid lines) and in the Hartle-Hawking state (dashed line). To appreciate the 
difference in F2 taken in the Hartle-Hawking and Unruh states, we show it in the inset. The G2 is the same for both states. 
The sonic horizon is located at zh = 0.025 cm. 



quadratic in the phonons field. To evaluate (pb,vb) up 
to 0(h) terms, the r.h.s. of the backreaction equations 
© and (fTUl) needs just to be evaluated on the classical 
background (p, v) of Sec. [H] 

The quantum state of the field in which the expectation 
values have to be computed depends on the physical situ- 
ation one wants to describe. For an isolated hole, the es- 
caping phonons radiation leads to a time variation of the 
underlying medium, i.e. pg(t,z) and Vs{t,z). The ap- 
propriate quantum state in this case is the analogue of the 
Unruh state |2j . In case the system is maintained in ther- 
mal equilibrium with the surroundings (that is, putting 
a sonic cavity in the subsonic asymptotic right region), 
the quantum state is the analogue of the Hartle-Hawking 
state the thermal equilibrium state at T = Th- In 
this case the system remains stationary, i.e. Pb(z) and 
v B (z). 

(2) 

Neglecting backscattering of the phonons, (T^J) can 



be approximated with the Polyakov stress tensor [9J, . 
Introducing for the sake of simplicity null coordinates 



x± 



dz 
c±v 



(11) 



the Polyakov stress tensor reads: 



and A±± are functions which depend on the choice of 
the quantum state of the phonons field. For the Unruh 
state: 



A++ = AY 



A__ 



, 

hk 2 

487TC 4 

For the Hartle-Hawking state instead 



(15) 
(16) 



A - aHH 



hk 2 

487TC 4 



(17) 



From Eqs. I|15|) - (|17[) it follows that, in the asymp- 
totic subsonic region z — > +00, A^ ± describes a flux 
of phonons at a temperature Th , whereas A™ describes 
a two-dimensional gas of phonons at thermal equilibrium 
at the temperature Th- To find first order in h correc- 
tions to the classical sonic black hole fluid configuration 
(p(z),v(z)) described in Eqs. JSJ and © we write 



PB 



p(z) - 



epi(t,z) 



(18) 
(19) 



where 



{TH) 



(T 



-12^ C C±± +A±± 



(2)\ _ 



f C-\\nC) 

D7T 



c 



pc* 



(12) 
(13) 

(14) 



with vb = d z ipB and e is a dimensionless expansion pa- 
rameter Ua: 



(20) 



\D\A H ' 

For our system e = 1.317 x 10 -14 . 



4 




0.03 0.032 



FIG. 4: Unruh state: time evolution of the backreaction equations for cut <C 1. Snapshots of the quantum correction to the 
velocity vi (left panel) and to the density (right panel) for an evolution time t on( j = 0.2t max ~ 2.09 x 10~ 2 ps (cut — 0.2). The 
vertical dashed line corresponds to the location of the classical sonic horizon. The delay between one snapshot and the other 
(between t = and t = i end ) is At ps 1.74 x 10~ 3 ps. 



The backreaction equations linearized in e then become 



IV. UNRUH STATE 



c 2 d z 



(2), 



(c — v) 2 (c + v) 2 



A + vvt H pi 

P 



— - — = £(j2 



eF 2 , (21) 
(22) 



Using the background equations QJ and ©, satisfied by 
p and v, the continuity equation can be rewritten as 



/ v v ' P v ' ,1 til F"2 
Pi + vp-y H —pi Vi +Wi = 7 

whereas the Bernoulli equation is 



ipx + «Vi 4 Pi = -r 



(23) 



(24) 



with a prime indicating derivative with respect to z. 

The profiles for the quantum sources Fi and Gi are 
depicted in Fig. [3] for the Unruh state (solid line) and for 
the Hartle-Hawking state (dashed line). The difference 
between the states is reflected on Fi only (and is shown in 
the inset in the left panel of the figure), while G2, being 
related to the trace anomaly which is state-independent, 
is unchanged. One can note the appearance of a maxi- 
mum and a minimum in the region z € [0.02, 0.03] cm. 
Outside this range, F2 and G2 rapidly drop to zero. The 
analysis of Ref. 6] , being limited to the region very near 
to zh , could not catch this non trivial structure. 



As said before, in Ref. @ the backreaction equations 
were analytically solved just for z « zj to allow a Taylor 
expansion of the sources up to linear terms. In this sec- 
tion we compute the numerical solution all over the noz- 
zle. We finite-difference the system of Eqs. (|23f) and l|24|) 
and solve it numerically in the time domain as an ini- 
tial value problem. The equations are discretized on an 
evenly spaced grid (0, z c ) with z c = 0.05 cm. Following a 
standard convention in numerical fluid mechanics |l3j| , we 
have used a staggered grid, i.e. both z — and z — z c are 
thought to lie on cell interfaces while the hydrodynam- 
ics quantities are defined on cell centers. As a result, the 
first point of our computational domain is z\ = Az/2 and 
the last is Zi max = z c — Az/2. We notice that Az is cho- 
sen so that the horizon is located at a cell interface. The 
reason for this is that, even though the square bracket on 
the r.h.s. of Eq. (|21|l is analytically regular at z = zh 
(where v — — c), the presence of the combination (c 4- v) 
at denominator in Eq. I|23|l can give problems (i.e., di- 
vision by zero) due to the discretization procedure. The 
use of a staggered grid bypasses this difficulty since the 
sonic point turns out to be always displaced with respect 
to the grid points. 

Initial conditions are chosen so that at t — the solu- 
tion is the classical one; i.e., p\(t = 0) = ipi(t = 0) = 0. 
Then the backreaction is switched on. As in Ref. 0, 
since the quantum sources are computed only for the 
static classical background, the validity of the solution 
is limited by the condition cut 1, where we have in- 
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FIG. 5: Profile of the quantum correction to the velocity v± (left panel) and to pi (right panel) due to backreaction in the 
Hartle-Hawking state (black lines), compared with the profile of vi and pi at t = i cn d in the Unruh state (blue lines). The 
vertical dashed lines correspond to the position of the sonic horizon at t — 0. 



troduced the constant k as 



dv 
dz 



(25) 



with dimension [length] _1 . For the sonic black hole con- 
sidered here, the short time condition determines a max- 
imum evolution time (cKt = 1) of t max = 0.104 [is; so 
it is possible to extract only informations about how the 
backreaction starts. 

Before discussing our numerical results, we briefly de- 
scribe the numerical algorithms implemented, further de- 
tails can be found in Appendix 1X1 

We are dealing with a system of Partial Differen- 
tial Equations (PDEs), where the equation for p\ is of 
convection-diffusion type, due to the parabolic term pro- 
portional to ip", while the equation for ip\ is a simple 
hyperbolic advection equation. As a result, the nu- 
merical algorithm must be designed accordingly |l2j . 
For the equation for ipx a simple first-order upwind 
method is well suited to solve it; for the parabolic 
equation we have implemented standard Forward-Timc- 
Centered-Space (FTCS) explicit method, as well as stan- 
dard Backward-Time-Centered-Space (BTCS) implicit 
method. Due to the short evolution time needed, the 
limitation of the time-step required by the FTCS and 
the consequent high number of iterations is not a draw- 
back; in any case, we tested one method versus the other 
and we obtained equivalent results. In fact, to have a sta- 
ble evolution, the time step is selected according to the 
condition At = aAz 2 /max(p), since p is the coefficient 
of the ip'{ term in the equation for pi. In addition, for 
the nozzle considered, we have checked that a resolution 
of Az — 2.5 x 10~ 5 cm (which corresponds to 2000 points 



covering the numerical domain) is sufficient to be in the 
convergence regime (see Appendix ^ for discussion) . 

In Fig. 2| we have snapshots of the time evolution of 
the profiles of v\ (left panel) and p\ (right panel). For 
this particular computation, we have considered a to- 
tal evolution time £ e nd = 0.2 £ max . The initial and fi- 
nal snapshots are depicted in red and blue respectively. 
The time delay between one snapshot and the other is 
At ~ 1.74 x 10~ 3 ps. The quantum corrected velocity is 
obtained as v\ = d z i\)\, the derivative being computed 
directly from the numerical data by means of a second 
order finite-difference approximation. 

The numerical solution confirms the near horizon be- 
haviour obtained in Ref. 0: the fluid slows down close 
to the horizon (v\ > 0, remember that v < because 
the fluid flows from right to left), causing the horizon to 
move to the left, and the total density decreases (pi < 0). 
In addition, now (even if for small times) it is possible 
to see the influence of the quantum corrections all over 
the sonic hole, and not just in the neighborhood of the 
horizon. As a consequence of the shape of the quantum 
sources F 2 and Gi (see Fig. |3J) the complex structure of 
Fig. ^emerges. One can see that in the region near the 
horizon the fluid slows down, but there are also regions 
where the phonons emission induces acceleration. 



V. HARTLE-HAWKING STATE 

The thermal equilibrium configuration of the Hartle- 
Hawking state is much simpler to treat. Since the time 
dependence drops off, the backreaction equations i|23|) 
and 124|) become a simple system of algebraic equations 
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relating p\ and v\. 

A (piv + pvi) = ( d£F 2 (0 + const 

Jzh 

G 2 . 



A vvi H pi 

P 



(26) 
(27) 



The integration constant in Eq. (|26|l is chosen to be zero 
in order to make the solution non singular on the horizon. 

The profile for v\ and p\ are depicted in Fig. CD (black 
line); for the sake of comparison, we show in the same 
plot the profile of v\ and p\ in the Unruh state for t = 
tend (blue line). In both cases the quantum backreaction 
correction to the velocity is positive at z = zh (vertical 
dashed line). 

In the region very close to the horizon one can make 
as in Ref. |(| a Taylor expansion for the background up 
to order 0((z — zh) 5 )- This allows the source terms to 
be evaluated up to linear term 



3 r 



F 2 = 



\D\A h k 

96tt 
tt 4 + 25tt 2 



- (tt 2 + 10) 

- 24 



(28) 



k{z - z H ) + 0{k\z - z H y) 



Go 



A 2 h c 2 k 2 
48tt 



x [(it 2 +6)k(z-z h ) + 0(k 2 (z-z h ) 2 )] 



(29) 



The corresponding quantum corrections to the velocity 
and to the density are 



AjjCK 2 
1927T _ 

52 + 35tt 2 



-k{z~z h )+0(k 2 (z-z h ) 2 ) 



(30) 



Pi 



192C7T _ 

44 - 19tt 2 



-k(z~z h )+0(k 2 (z-z h ) 2 ) 



(31) 

Setting vb = V + ev% = —c one finds the quantum 
corrected position of the horizon z q H 



-H 



ZH 



1927T 



eA H n , 



(32) 



which is shifted to the left of zh- 

The quantum corrected equilibrium temperature can also 



be simply obtained by evaluating Eq. JSJ at 
v replaced by Vb- The result is 



-H 



ml 
1 H 



flCK 
2,7TK B 



1 



cA 



H 



768tt 



52 



35tt 2 



4\ 2 
7T K 



with 



(33) 



which indicates that, taking into account the backreac- 
tion, the equilibrium temperature is lowered. 



VI. CONCLUSIONS 

Using the continuity and Bernoulli equations, the 
quantum correction (first order in h) to a classical sta- 
tionary flow describing an acoustic black hole has been 
evaluated in a one-dimensional approximation. 
The quantum corrections to the velocity V\ and to 
the density p± profiles for the equilibrium configura- 
tion (Hartle-Hawking state) are depicted in Fig. [S] (black 
lines) . The phonons backreaction causes the fluid to slow 
down in the supersonic region, with the consequence of 
a shift of the horizon to the left of the waist of the noz- 
zle (see Eq. (|32p'l. In the subsonic region the velocity 
increases, but the magnitude of the change is smaller 
than the previous one. One finds a similar shape for the 
density correction p±, which increases in the supersonic 
region and slightly decreases in the subsonic one. Fi- 
nally the equilibrium temperature appears to have been 
lowered by the backreaction from its zero-order value 
Hck/2ttkb (see Eq. 

For the time-dependent case (Unruh state) the analysis 
has been restricted to very short times after switching on 
the phonons radiation. This because the quantum source 
(F 2 and G 2 in Eqs. I|21() - (|22|l ') has been computed only for 
the classical background, which just represents the initial 
configuration of the acoustic black hole. A more rigorous 
analysis requires the time-dependence of the sources to 
be included. 

Within these limitations, one sees (Fig. |2J a decelera- 
tion of the fluid in the supersonic region, which causes a 
drift of the horizon towards the left of the nozzle. Two 
acceleration regions also appear on both sides of the hori- 
zon, but the intensity of the effect is lower. On the other 
hand, the density correction p\ reflects the behaviour it 
shows in the Hartle-Hawking state (on a reduced scale). 
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APPENDIX A: NUMERICAL SCHEMES 

In this section we report explicitly the time-evolution 
algorithms. In the main text, we said that we used a 
standard upwind method for the equation for tpi and 
standard explicit Forward-Time-Centered-Space (FTCS) 
or an implicit Backward-Time-Centered-Space (BTCS) 
schemes for that for p±. In practice the upwind method 
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reads [3 The FTCS (explicit) and the BTCS (implicit) schemes 

for the equation for p\ respectively read 

r^t =^-^(^+1-^,0 

+ ^t(~pl i + G 3>i A^A . (Al) 



ViAt 
'~2Az 



n+l 
PlA 



n+1 
PlA-1 



; Pl,i 

At 
1 



2A 



z 

2 / 



2 / 



At p"; 



ViAt 



Tt + l 



2A; 



~ Pi. 



Az 2 



2Az 



Az 2 



(A2) 



(A3) 



where the cell index i runs from one to i max - In the case 
of the BTCS scheme p\ is obtained at every time slice 
(labelled by index n) as the solution of a tridiagonal lin- 
ear system of the form + + CiU™ +1 = /" that 
can be accomplished by a standard lower-upper (LU) de- 
composition of the matrix to be inverted 14]. A careful 
treatment of the boundaries 2 = and z = z c of the 
numerical domain is crucial for selecting the correct so- 
lution, especially when the implicit method is employed 
and so the inversion of the associated coefficient matrix is 
concerned. According to the physical meaning of the Un- 
ruh state, we impose outgoing conditions at both bound- 
aries, i.e. = it™ at % = 1 and = u™ at i = i max 
where u n can be either p\ or tpi . If p\ is solved using the 
FTCS scheme, since the method is explicit and no matrix 
inversion is needed, the problem of setting correct bound- 
ary conditions is less important; in fact, it is enough to 
put the boundaries far enough from zh to avoid any in- 
fluence on the evolution. For this kind of equations, sta- 
bility has also proved to be an issue. Implementing the 
FTCS scheme, to have a stable evolution the time step is 



selected according to the condition At = aAz 2 /max(p). 
For the nozzle model discussed in this paper, we have 
used a = 1 X 10 -4 to avoid stability problems and the 
same choice was kept also for the BTCS scheme, which re- 
sults in roughly 3 x 10 4 integration steps. This is not par- 
ticularly expensive from the computational point of view. 
For example, for the results presented here we used a res- 
olution of Az = 2.5 x 10~ 5 cm (corresponding to 2000 
grid points) and the total evolution took ~ lis time evo- 
lution on a single-processor machine with a Pentium™ 
M processor at 1.3GHz. The code was compiled using an 
Intel™ Fortran Compiler. 

We checked convergence of the numerical method (us- 
ing both BTCS or FTCS schemes) using resolutions of 
500, 1000, 2000 and 4000 points, considering the case of 
8000 points as reference. We computed the error A/ with 
respect to the reference resolution as a root mean square 
for / = tpi and f = p\. From the relation A/ = K.Az a we 
evaluated the convergence rate a and we obtained a ~ 1.3 
for both Vi and p x . We have verified that 2000 4000 
grid points are sufficient to be in the convergence regime. 
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